/*
 * =====================================================================================
 *
 *       Filename:  loopybp.h
 *
 *    Description:  loopy belief propagation
 *
 *        Version:  1.0
 *        Created:  2009年06月29日 10时34分06秒
 *       Revision:  none
 *       Compiler:  gcc
 *
 *         Author:  Ying Wang (WY), ywang@nlpr.ia.ac.cn
 *        Company:  Institute of Automation, Chinese Academy of Sciences
 *
 * =====================================================================================
 */

#include "ncmatrix.h"
#include <vector>
#define MAXD 1e20
#define POTTS 1000
class preMRF
{
public:
    class Mat_mm : public NCmatrix<double>
    {
    public:
        Mat_mm()
            :NCmatrix<double> (mmwidth,mmheight)
        {
        }

    };
    preMRF(int m,int n)
    {
        mmwidth = m;
	mmheight = n;
    }
    static int mmwidth;
    static int mmheight;

};
int preMRF::mmwidth  = -1;
int preMRF::mmheight = -1;

class MRF : public preMRF
{
public:
	MRF(int row, int column, int mlabels,int iter=8)
		: preMRF(row,column)
		, m_row(row)
		, m_column(column)
		, m_labels(mlabels)
		, iter_num(iter)
		, leftm(m_labels)
		, rightm(m_labels)
		, upm(m_labels)
		, downm(m_labels)
		, label(m_row,m_column)
	{
	}
	void maxprodbp();
	void fastmaxprodbp();
	
public:
	//void computeLeftMsg();
	virtual double dataitem(int pi,int pj,int labelp) const = 0;
	virtual double prioritem(int pi,int pj,int qi,int qj, int labelp,int labelq) const =0;
	void initialMsg();
public:
	inline NCmatrix<int> labelmatrix()const
	{
		return label;
	}
private:
	int m_row;
	int m_column;
	int m_labels;
	int iter_num;
	std::vector<Mat_mm> leftm; // left message 
	std::vector<Mat_mm> rightm; // right message
	std::vector<Mat_mm> upm; // up message
	std::vector<Mat_mm> downm; // down message
	NCmatrix<int> label;
};

void MRF::fastmaxprodbp()
{
	initialMsg();
	std::cout<<"fast loopy min-sum algorithm"<<std::endl;
	int c,i,j,z;
	double enleft,enup,enright,endown,leftmin,upmin,rightmin,downmin;
	std::vector<Mat_mm> last_leftm,last_rightm,last_upm,last_downm;
	for( int iter=0;iter<iter_num; iter++ )
	{
		last_leftm  = leftm;
		last_rightm = rightm;
		last_upm    = upm;
		last_downm  = downm;
		std::cout<<"iter number: " << iter+1<<std::endl;
		//
		i=0; j=0;
		rightmin = MAXD;
		downmin = MAXD;
		for( c=0; c<m_labels; c++ )
		{
			enright = dataitem(i,j+1,c) + last_upm[c][i][j+1] + last_rightm[c][i][j+1]
				+ last_downm[c][i][j+1];
			endown  = dataitem(i+1,j,c) + last_leftm[c][i+1][j] + last_rightm[c][i+1][j]
				+ last_downm[c][i+1][j];
			if(enright<rightmin)
				rightmin = enright;
			if(endown<downmin)
				downmin = endown;
		}
		for( z=0; z<m_labels; z++ )
		{
			enright = dataitem(i,j+1,z) + last_upm[z][i][j+1] + last_rightm[z][i][j+1]
				+ last_downm[z][i][j+1];
			endown  = dataitem(i+1,j,z) + last_leftm[z][i+1][j] + last_rightm[z][i+1][j]
				+ last_downm[z][i+1][j];
			leftm[z][i][j] = 0;
			rightm[z][i][j] = (enright<(POTTS+rightmin) ? enright : (POTTS+rightmin));
			upm[z][i][j]    = 0;
			downm[z][i][j]  = (endown<(POTTS+downmin) ? endown : (POTTS+downmin));	
			
		}
		//
		i=0; j = m_column-1;
		leftmin = MAXD;
		downmin = MAXD;
		for( c=0; c<m_labels; c++ )
		{
			enleft = dataitem(i,j-1,c) + last_leftm[c][i][j-1] + last_upm[c][i][j-1]
				+ last_downm[c][i][j-1];
			endown  = dataitem(i+1,j,c) + last_leftm[c][i+1][j] + last_rightm[c][i+1][j]
				+ last_downm[c][i+1][j];
			if(enleft<leftmin)
				leftmin = enleft;
			if(endown<downmin)
				downmin = endown;
		}
		for( z=0; z<m_labels; z++ )
		{
			enleft = dataitem(i,j-1,z) + last_leftm[z][i][j-1] + last_leftm[z][i][j-1] 
				+ last_downm[z][i][j-1];

			endown  = dataitem(i+1,j,z) + last_leftm[z][i+1][j] + last_rightm[z][i+1][j]
				+ last_downm[z][i+1][j];

			leftm[z][i][j]  = (enleft<(POTTS+leftmin) ? enleft : (POTTS+leftmin));
			rightm[z][i][j] = 0;
			upm[z][i][j]    = 0;
			downm[z][i][j]  = (endown<(POTTS+downmin) ? endown : (POTTS+downmin));	
			
		}

		//
		i=m_row-1; j=0;
		rightmin = MAXD;
		upmin = MAXD;
		for( c=0; c<m_labels; c++ )
		{

			enup   = dataitem(i-1,j,c) + last_leftm[c][i-1][j] + last_upm[c][i-1][j]
				+ last_rightm[c][i-1][j];
			enright = dataitem(i,j+1,c) + last_upm[c][i][j+1] + last_rightm[c][i][j+1]
				+ last_downm[c][i][j+1];

			if(enup<upmin)
				upmin = enup;
			if(enright<rightmin)
				rightmin = enright;
		}
		for( z=0; z<m_labels; z++ )
		{

			enup   = dataitem(i-1,j,z) + last_leftm[z][i-1][j] + last_upm[z][i-1][j]
				+ last_rightm[z][i-1][j];
			enright = dataitem(i,j+1,z) + last_upm[z][i][j+1] + last_rightm[z][i][j+1]
				+ last_downm[z][i][j+1];

			leftm[z][i][j] =0;
			rightm[z][i][j] = (enright<(POTTS+rightmin) ? enright : (POTTS+rightmin));
			upm[z][i][j]    = (enup<(POTTS+upmin) ? enup : (POTTS+upmin));
			downm[z][i][j] = 0;
			
		}

		//
		i=m_row-1; j=m_column-1;

		leftmin = MAXD;
		upmin = MAXD;
		for( c=0; c<m_labels; c++ )
		{
			enleft = dataitem(i,j-1,c) + last_leftm[c][i][j-1] + last_upm[c][i][j-1]
				+ last_downm[c][i][j-1];
			enup   = dataitem(i-1,j,c) + last_leftm[c][i-1][j] + last_upm[c][i-1][j]
				+ last_rightm[c][i-1][j];
			if(enleft<leftmin)
				leftmin = enleft;
			if(enup<upmin)
				upmin = enup;
		}
		for( z=0; z<m_labels; z++ )
		{
			enleft = dataitem(i,j-1,z) + last_leftm[z][i][j-1] + last_leftm[z][i][j-1] 
				+ last_downm[z][i][j-1];
			enup   = dataitem(i-1,j,z) + last_leftm[z][i-1][j] + last_upm[z][i-1][j]
				+ last_rightm[z][i-1][j];

			leftm[z][i][j]  = (enleft<(POTTS+leftmin) ? enleft : (POTTS+leftmin));
			rightm[z][i][j] =0;
			upm[z][i][j]    = (enup<(POTTS+upmin) ? enup : (POTTS+upmin));
			downm[z][i][j] = 0;
		}
		//first row
		i=0;
		for( j=1; j<m_column-1; j++ )
		{
			leftmin = MAXD;
			//upmin   = MAXD;
			rightmin= MAXD;
			downmin = MAXD;
			for( c=0; c<m_labels; c++ )
			{
				enleft = dataitem(i,j-1,c) + last_leftm[c][i][j-1] + last_upm[c][i][j-1]
					+ last_downm[c][i][j-1];
				//enup   = dataitem(i-1,j,c) + last_leftm[c][i-1][j] + last_upm[c][i-1][j]
				//	+ last_rightm[c][i-1][j];
				enright = dataitem(i,j+1,c) + last_upm[c][i][j+1] + last_rightm[c][i][j+1]
					+ last_downm[c][i][j+1];
				endown  = dataitem(i+1,j,c) + last_leftm[c][i+1][j] + last_rightm[c][i+1][j]
					+ last_downm[c][i+1][j];
				if(enleft<leftmin)
					leftmin = enleft;
				//if(enup<upmin)
				//	upmin = enup;
				if(enright<rightmin)
					rightmin = enright;
				if(endown<downmin)
					downmin = endown;
			}
			for( z=0; z<m_labels; z++ )
			{
				enleft = dataitem(i,j-1,z) + last_leftm[z][i][j-1] + last_leftm[z][i][j-1] 
					+ last_downm[z][i][j-1];
				//enup   = dataitem(i-1,j,z) + last_leftm[z][i-1][j] + last_upm[z][i-1][j]
				//	+ last_rightm[z][i-1][j];
				enright = dataitem(i,j+1,z) + last_upm[z][i][j+1] + last_rightm[z][i][j+1]
					+ last_downm[z][i][j+1];
				endown  = dataitem(i+1,j,z) + last_leftm[z][i+1][j] + last_rightm[z][i+1][j]
					+ last_downm[z][i+1][j];

				leftm[z][i][j]  = (enleft<(POTTS+leftmin) ? enleft : (POTTS+leftmin));
				rightm[z][i][j] = (enright<(POTTS+rightmin) ? enright : (POTTS+rightmin));
				upm[z][i][j]    = 0;
				downm[z][i][j]  = (endown<(POTTS+downmin) ? endown : (POTTS+downmin));	
				
			}
		}
		//last_row
		i = m_row-1;
		for( j=1; j<m_column-1; j++ )
		{
			leftmin = MAXD;
			upmin   = MAXD;
			rightmin= MAXD;
			//downmin = MAXD;
			for( c=0; c<m_labels; c++ )
			{
				enleft = dataitem(i,j-1,c) + last_leftm[c][i][j-1] + last_upm[c][i][j-1]
					+ last_downm[c][i][j-1];
				enup   = dataitem(i-1,j,c) + last_leftm[c][i-1][j] + last_upm[c][i-1][j]
					+ last_rightm[c][i-1][j];
				enright = dataitem(i,j+1,c) + last_upm[c][i][j+1] + last_rightm[c][i][j+1]
					+ last_downm[c][i][j+1];
			//	endown  = dataitem(i+1,j,c) + last_leftm[c][i+1][j] + last_rightm[c][i+1][j]
			//		+ last_downm[c][i+1][j];
				if(enleft<leftmin)
					leftmin = enleft;
				if(enup<upmin)
					upmin = enup;
				if(enright<rightmin)
					rightmin = enright;
			//	if(endown<downmin)
			//		downmin = endown;
			}
			for( z=0; z<m_labels; z++ )
			{
				enleft = dataitem(i,j-1,z) + last_leftm[z][i][j-1] + last_leftm[z][i][j-1] 
					+ last_downm[z][i][j-1];
				enup   = dataitem(i-1,j,z) + last_leftm[z][i-1][j] + last_upm[z][i-1][j]
					+ last_rightm[z][i-1][j];
				enright = dataitem(i,j+1,z) + last_upm[z][i][j+1] + last_rightm[z][i][j+1]
					+ last_downm[z][i][j+1];
		//		endown  = dataitem(i+1,j,z) + last_leftm[z][i+1][j] + last_rightm[z][i+1][j]
		//			+ last_downm[z][i+1][j];

				leftm[z][i][j]  = (enleft<(POTTS+leftmin) ? enleft : (POTTS+leftmin));
				rightm[z][i][j] = (enright<(POTTS+rightmin) ? enright : (POTTS+rightmin));
				upm[z][i][j]    = (enup<(POTTS+upmin) ? enup : (POTTS+upmin));
				downm[z][i][j]   = 0;
		//		downm[z][i][j]  = (endown<(POTTS+downmin) ? endown : (POTTS+downmin));	
				
			}
		}
		
		//first column
		j = 0;
		for( i=1; i<m_row-1; i++ )
		{
			//leftmin = MAXD;
			upmin   = MAXD;
			rightmin= MAXD;
			downmin = MAXD;
			for( c=0; c<m_labels; c++ )
			{
				//enleft = dataitem(i,j-1,c) + last_leftm[c][i][j-1] + last_upm[c][i][j-1]
				//	+ last_downm[c][i][j-1];
				enup   = dataitem(i-1,j,c) + last_leftm[c][i-1][j] + last_upm[c][i-1][j]
					+ last_rightm[c][i-1][j];
				enright = dataitem(i,j+1,c) + last_upm[c][i][j+1] + last_rightm[c][i][j+1]
					+ last_downm[c][i][j+1];
				endown  = dataitem(i+1,j,c) + last_leftm[c][i+1][j] + last_rightm[c][i+1][j]
					+ last_downm[c][i+1][j];
			//	if(enleft<leftmin)
			//		leftmin = enleft;
				if(enup<upmin)
					upmin = enup;
				if(enright<rightmin)
					rightmin = enright;
				if(endown<downmin)
					downmin = endown;
			}
			for( z=0; z<m_labels; z++ )
			{
			//	enleft = dataitem(i,j-1,z) + last_leftm[z][i][j-1] + last_leftm[z][i][j-1] 
			//		+ last_downm[z][i][j-1];
				enup   = dataitem(i-1,j,z) + last_leftm[z][i-1][j] + last_upm[z][i-1][j]
					+ last_rightm[z][i-1][j];
				enright = dataitem(i,j+1,z) + last_upm[z][i][j+1] + last_rightm[z][i][j+1]
					+ last_downm[z][i][j+1];
				endown  = dataitem(i+1,j,z) + last_leftm[z][i+1][j] + last_rightm[z][i+1][j]
					+ last_downm[z][i+1][j];

				leftm[z][i][j] = 0;
		//		leftm[z][i][j]  = (enleft<(POTTS+leftmin) ? enleft : (POTTS+leftmin));
				rightm[z][i][j] = (enright<(POTTS+rightmin) ? enright : (POTTS+rightmin));
				upm[z][i][j]    = (enup<(POTTS+upmin) ? enup : (POTTS+upmin));
				downm[z][i][j]  = (endown<(POTTS+downmin) ? endown : (POTTS+downmin));	
				
			}
		}

		//last column

		j = m_column-1;
		for( i=1; i<m_row-1; i++ )
		{
			leftmin = MAXD;
			upmin   = MAXD;
			//rightmin= MAXD;
			downmin = MAXD;
			for( c=0; c<m_labels; c++ )
			{
				enleft = dataitem(i,j-1,c) + last_leftm[c][i][j-1] + last_upm[c][i][j-1]
					+ last_downm[c][i][j-1];
				enup   = dataitem(i-1,j,c) + last_leftm[c][i-1][j] + last_upm[c][i-1][j]
					+ last_rightm[c][i-1][j];
				//enright = dataitem(i,j+1,c) + last_upm[c][i][j+1] + last_rightm[c][i][j+1]
				//	+ last_downm[c][i][j+1];
				endown  = dataitem(i+1,j,c) + last_leftm[c][i+1][j] + last_rightm[c][i+1][j]
					+ last_downm[c][i+1][j];
				if(enleft<leftmin)
					leftmin = enleft;
				if(enup<upmin)
					upmin = enup;
				//if(enright<rightmin)
				//	rightmin = enright;
				if(endown<downmin)
					downmin = endown;
			}
			for( z=0; z<m_labels; z++ )
			{
				enleft = dataitem(i,j-1,z) + last_leftm[z][i][j-1] + last_leftm[z][i][j-1] 
					+ last_downm[z][i][j-1];
				enup   = dataitem(i-1,j,z) + last_leftm[z][i-1][j] + last_upm[z][i-1][j]
					+ last_rightm[z][i-1][j];
				//enright = dataitem(i,j+1,z) + last_upm[z][i][j+1] + last_rightm[z][i][j+1]
				//	+ last_downm[z][i][j+1];
				endown  = dataitem(i+1,j,z) + last_leftm[z][i+1][j] + last_rightm[z][i+1][j]
					+ last_downm[z][i+1][j];

				leftm[z][i][j]  = (enleft<(POTTS+leftmin) ? enleft : (POTTS+leftmin));
				rightm[z][i][j] = 0.;
				//rightm[z][i][j] = (enright<(POTTS+rightmin) ? enright : (POTTS+rightmin));
				upm[z][i][j]    = (enup<(POTTS+upmin) ? enup : (POTTS+upmin));
				downm[z][i][j]  = (endown<(POTTS+downmin) ? endown : (POTTS+downmin));	
				
			}
		}

		//M[1:m_row-2]
		for( i=1; i<m_row-1; i++ )
		{
			for( j=1; j<m_column-1; j++ )
			{
				leftmin = MAXD;
				upmin   = MAXD;
				rightmin= MAXD;
				downmin = MAXD;
				for( c=0; c<m_labels; c++ )
				{
					enleft = dataitem(i,j-1,c) + last_leftm[c][i][j-1] + last_upm[c][i][j-1]
						+ last_downm[c][i][j-1];
					enup   = dataitem(i-1,j,c) + last_leftm[c][i-1][j] + last_upm[c][i-1][j]
						+ last_rightm[c][i-1][j];
					enright = dataitem(i,j+1,c) + last_upm[c][i][j+1] + last_rightm[c][i][j+1]
						+ last_downm[c][i][j+1];
					endown  = dataitem(i+1,j,c) + last_leftm[c][i+1][j] + last_rightm[c][i+1][j]
						+ last_downm[c][i+1][j];
					if(enleft<leftmin)
						leftmin = enleft;
					if(enup<upmin)
						upmin = enup;
					if(enright<rightmin)
						rightmin = enright;
					if(endown<downmin)
						downmin = endown;
				}
				for( z=0; z<m_labels; z++ )
				{
					enleft = dataitem(i,j-1,z) + last_leftm[z][i][j-1] + last_leftm[z][i][j-1] 
						+ last_downm[z][i][j-1];
					enup   = dataitem(i-1,j,z) + last_leftm[z][i-1][j] + last_upm[z][i-1][j]
						+ last_rightm[z][i-1][j];
					enright = dataitem(i,j+1,z) + last_upm[z][i][j+1] + last_rightm[z][i][j+1]
						+ last_downm[z][i][j+1];
					endown  = dataitem(i+1,j,z) + last_leftm[z][i+1][j] + last_rightm[z][i+1][j]
						+ last_downm[z][i+1][j];

					leftm[z][i][j]  = (enleft<(POTTS+leftmin) ? enleft : (POTTS+leftmin));
					rightm[z][i][j] = (enright<(POTTS+rightmin) ? enright : (POTTS+rightmin));
					upm[z][i][j]    = (enup<(POTTS+upmin) ? enup : (POTTS+upmin));
					downm[z][i][j]  = (endown<(POTTS+downmin) ? endown : (POTTS+downmin));	

				}
			}
		}
	}
	//
	double tmpmin,tmp;
	int zmin;
	for( i=0; i<m_row; i++ )
	{
		for (j=0; j<m_column; j++ )
		{
			tmpmin = 1e20;
			for( z=0; z<m_labels; z++ )
			{
				tmp = dataitem(i,j,z) + leftm[z][i][j] + 
					rightm[z][i][j] + upm[z][i][j] + downm[z][i][j];
				if(tmp<tmpmin)
				{
					tmpmin = tmp;
					zmin = z;
				}
			}
			label[i][j] = zmin;
		}
	}
}
void MRF::maxprodbp()
{
	initialMsg();
	std::cout<<"loopy min sum algorithm"<<std::endl;
	int c,i,j,z;
	double enleft,enup,enright,endown,leftmin,upmin,rightmin,downmin;
	std::vector<Mat_mm> last_leftm, last_rightm, last_upm, last_downm ;
	for( int iter=0;iter<iter_num;iter++)
	{
		last_leftm = leftm;
		last_rightm = rightm;
		last_upm = upm;
		last_downm = downm;
		std::cout<<"iter number: "<<iter+1<<std::endl;

		for( i=0; i<m_row; i++ )
		{
			for( j=0; j<m_column; j++ )
			{
				//std::cout<<"(i,j)= "<<i<<","<<j<<std::endl;
				for( z=0; z<m_labels; z++ )
				{
					if(i==0 && j==0 )
					{
						upm[z][i][j] = 0.;
						leftm[z][i][j] = 0.;
					}
					if (i==0 && j== m_column-1)
					{
						upm[z][i][j] = 0.;
						rightm[z][i][j] =0.;
					}
					if (i==m_row-1 && j== 0 )
					{
						leftm[z][i][j] =0.;
						downm[z][i][j] = 0.;
					}
					if (i==m_row-1 && j== m_column-1)
					{
						rightm[z][i][j] =0.;
						downm[z][i][j] = 0.;
					}
					if ( i==0 && j>0 && j<m_column-1 )
					{
						upm[z][i][j] = 0.;
						leftmin = MAXD;
						rightmin = MAXD;
						downmin = MAXD;
						for( c=0; c<m_labels; c++ )
						{
							enleft = dataitem(i,j-1,c) 
								 + prioritem(i,j-1,i,j,c,z)
								 + last_leftm[c][i][j-1]
								 + last_downm[c][i][j-1];
							enright = dataitem(i,j+1,c)
								  + prioritem(i,j+1,i,j,c,z)
								  + last_rightm[c][i][j+1]
								  + last_downm[c][i][j+1];
							endown = dataitem(i+1,j,c)
								 + prioritem(i+1,j,i,j,c,z)
								 + last_leftm[c][i+1][j]
								 + last_rightm[c][i+1][j]
								 + last_downm[c][i+1][j];
							if(enleft<leftmin)
								leftmin = enleft;
							if(enright<rightmin)
								rightmin = enright;
							if(endown<downmin)
								downmin = endown;
								 
						}
						leftm[z][i][j] = leftmin;
						rightm[z][i][j] = rightmin;
						downm[z][i][j] = downmin;

					}
					if( i==m_row-1 && j>0 && j<m_column-1 )
					{
						downm[z][i][j] = 0.;
						leftmin = MAXD;
						rightmin = MAXD;
						upmin = MAXD;
						for( c=0; c<m_labels; c++ )
						{
							enleft = dataitem(i,j-1,c) 
								 + prioritem(i,j-1,i,j,c,z)
								 + last_leftm[c][i][j-1]
								 + last_downm[c][i][j-1];
							enright = dataitem(i,j+1,c)
								  + prioritem(i,j+1,i,j,c,z)
								  + last_rightm[c][i][j+1]
								  + last_downm[c][i][j+1];
							enup = dataitem(i-1,j,c)
								 + prioritem(i-1,j,i,j,c,z)
								 + last_leftm[c][i-1][j]
								 + last_rightm[c][i-1][j]
								 + last_upm[c][i-1][j];
							if(enleft<leftmin)
								leftmin = enleft;
							if(enright<rightmin)
								rightmin = enright;
							if(enup<upmin)
								upmin = enup;
								 
						}
						leftm[z][i][j] = leftmin;
						rightm[z][i][j] = rightmin;
						upm[z][i][j] = upmin;

					}
					if( j==0 && i>0 && i<m_row-1 )
					{
						leftm[z][i][j] = 0.;
						upmin = MAXD;
						rightmin = MAXD;
						downmin = MAXD;
						for( c=0; c<m_labels; c++ )
						{
							enup = dataitem(i-1,j,c) 
								 + prioritem(i-1,j,i,j,c,z)
								 + last_upm[c][i-1][j]
								 + last_rightm[c][i-1][j];
							enright = dataitem(i,j+1,c)
								  + prioritem(i,j+1,i,j,c,z)
								  + last_upm[c][i][j+1]
								  + last_rightm[c][i][j+1]
								  + last_downm[c][i][j+1];
							endown = dataitem(i+1,j,c)
								 + prioritem(i+1,j,i,j,c,z)
								 + last_rightm[c][i+1][j]
								 + last_downm[c][i+1][j];
							if(enup<upmin)
								upmin = enup;
							if(enright<rightmin)
								rightmin = enright;
							if(endown<downmin)
								downmin = endown;
								 
						}
						upm[z][i][j] = upmin;
						rightm[z][i][j] = rightmin;
						downm[z][i][j] = downmin;

					}
					if( j==m_column-1 && i>0 && i<m_row-1 )
					{
						rightm[z][i][j] = 0.;
						upmin = MAXD;
						leftmin = MAXD;
						downmin = MAXD;
						for( c=0; c<m_labels; c++ )
						{
							enup = dataitem(i-1,j,c) 
								 + prioritem(i-1,j,i,j,c,z)
								 + last_upm[c][i-1][j]
								 + last_leftm[c][i-1][j];
							enleft = dataitem(i,j-1,c)
								  + prioritem(i,j-1,i,j,c,z)
								  + last_upm[c][i][j-1]
								  + last_leftm[c][i][j-1]
								  + last_downm[c][i][j-1];
							endown = dataitem(i+1,j,c)
								 + prioritem(i+1,j,i,j,c,z)
								 + last_rightm[c][i+1][j]
								 + last_downm[c][i+1][j];
							if(enup<upmin)
								upmin = enup;
							if(enleft<leftmin)
								leftmin = enleft;
							if(endown<downmin)
								downmin = endown;
								 
						}
						upm[z][i][j] = upmin;
						leftm[z][i][j] = leftmin;
						downm[z][i][j] = downmin;

					}
					if( i>0 && i<m_row-1 && j>0 && j<m_column-1 )
					{
						leftmin  = MAXD;
						upmin    = MAXD;
						rightmin = MAXD;
						downmin  = MAXD;
						for( c=0; c<m_labels; c++ )
						{
							enleft  = dataitem(i,j-1,c)
								  + prioritem(i,j-1,i,j,c,z)
								  + last_leftm[c][i][j-1]
								  + last_upm[c][i][j-1]
								  + last_downm[c][i][j-1];
							enright = dataitem(i,j+1,c)
								  + prioritem(i,j+1,i,j,c,z)
								  + last_rightm[c][i][j+1]
								  + last_upm[c][i][j+1]
								  + last_downm[c][i][j+1];
							enup    = dataitem(i-1,j,c)
								  + prioritem(i-1,j,i,j,c,z)
								  + last_leftm[c][i-1][j]
								  + last_rightm[c][i-1][j]
								  + last_downm[c][i-1][j];
							endown  = dataitem(i+1,j,c)
								  + prioritem(i+1,j,i,j,c,z)
								  + last_leftm[c][i+1][j]
								  + last_rightm[c][i+1][j]
								  + last_downm[c][i+1][j];
							if(enup<upmin)
								upmin = enup;
							if(enleft<leftmin)
								leftmin = enleft;
							if(enright<rightmin)
								rightmin = enright;
							if(endown<downmin)
								downmin = endown;
						}
						upm[z][i][j] = upmin;
						leftm[z][i][j] = leftmin;
						rightm[z][i][j] = rightmin;
						downm[z][i][j] = downmin;
					}


				}
			}
		}
	}


	double tmpmin,tmp;
	int zmin;
	for( i=0; i<m_row; i++ )
	{
		for (j=0; j<m_column; j++ )
		{
			tmpmin = 1e20;
			for( z=0; z<m_labels; z++ )
			{
				tmp = dataitem(i,j,z) + leftm[z][i][j] + 
					rightm[z][i][j] + upm[z][i][j] + downm[z][i][j];
				if(tmp<tmpmin)
				{
					tmpmin = tmp;
					zmin = z;
				}
			}
			label[i][j] = zmin;
		}
	}
}

/*void MRF::computeLeftMsg()*/
/*{*/

/*}*/
void MRF::initialMsg()
{
int zi,i,j;
for(zi=0;zi<m_labels;zi++)
{
	for( i=0; i<m_row; i++ )
	{
		for( j=0; j<m_column; j++ )
		{
			upm[zi][i][j] = 0.;
			downm[zi][i][j] =0.;
			leftm[zi][i][j] = 0.;
			rightm[zi][i][j] = 0.;
		}
	}
	
}
}

/*inline double MRF::dataitem(int pi,int pj,int labelp)
{
	double dataenergy;
	return dataenergy;
}
inline double MRF::prioritem(int pi,int pj, int qi, int qj,int labelp,int labelq)
{
	double priorenergy;
	return priorenergy;
}*/
